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Generation of Porous Particle Structures using the Void 
Expansion Method 

Ov 

| Iwan Schenker • Frank T. Filser • Hans J. Herrmann • 
CNJ ■ Ludwig J. Gauckler 

> : 
o ■ 



o 

c3 



I 

c 

o 
o 



o 



o 



X 



Received: 29. 6. 2008 / Accepted: 26. 1. 2009 

Abstract The newly developed "void expansion 
method" allows for an efficient generation of porous 
packings of spherical particles over a wide range of 
volume fractions using the discrete element method. 
Particles are randomly placed under addition of much 
smaller " void-particles" . Then, the void-particle radius 
is increased repeatedly, thereby rearranging the struc- 
tural particles until formation of a dense particle pack- 
ing. 

The structural particles' mean coordination number 
was used to characterize the evolving microstructures. 
At some void radius, a transition from an initially low 
to a higher mean coordination number is found, which 

' was used to characterize the influence of the various 
simulation parameters. For structural and void-particle 
stiffnesses of the same order of magnitude, the transi- 
tion is found at constant total volume fraction slightly 
below the random close packing limit. For decreasing 
void-particle stiffness the transition is shifted towards 

, a smaller void-particle radius and becomes smoother. 

Keywords colloid • coordination number • discrete 
element method • microstructure generation • porosity 



1 Introduction 

Mechanical tests on coagulated colloids have shown 
that the local arrangement of the colloidal particles 
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has a strong influence on the macroscopic mechani- 
cal properties. Colloids with a more "heterogeneous" 
microstructure possess up to one order of magnitude 
higher elastic modulii and yield strengths than their 
"homogeneous" counterparts pQ. 

Experimentally, colloidal microstructures with dif- 
ferent degrees of heterogeneity are obtained using an 
internal gelation method (DCC = Direct Coagulation 
Casting [2J[3]). The method allows for an in-situ, i.e. 
undisturbed, transition of the inter-particle potential 
from repulsive to attractive. There are two principal 
pathways leading to different microstructures: chang- 
ing the pH of the suspension (ApH-method) or increas- 
ing the ionic strength in the suspension (AI- method). 
The first pathway shifts the pH to the particles' iso- 
electric point and produces more "homogeneous" mi- 
crostructures through diffusion limited aggregation. In 
the second pathway the ionic strength in the suspen- 
sion is increased at a constant pH which compresses 
the Debye length of the repulsive potential leading to 
more "heterogeneous" microstructures due to reaction 
rate limited aggregation of the particles [4] . 

Alternatively, heterogeneous microstructures can as 
well be obtained by ApH-destabilization in conjunction 
with small amounts of alkali-swellable polymer particles 
(ASP), 80 nm in diameter in the unswollen state [5]. 
The ASP particles were admixed to the structural par- 
ticles of 200 nm in diameter under acidic conditions and 
swelled upon increasing pH during the internal gelling 
reaction of the DCC process unfolding to 800 nm in 
diameter, and thus pushing the structural particles in 
their vicinity. Thereby, larger pores and thus more het- 
erogeneous microstructures are produced. Those sam- 
ples with ASP exhibit much higher mechanical proper- 
ties than samples without ASP. In particular, ASP sam- 
ples present comparably high mechanical properties as 
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samples with heterogeneous microstructures produced 
by the Al-method. 

These experimental findings suggest that the col- 
loid's microstructure strongly determines its macro- 
scopic mechanical properties. The relation between 
structure and mechanical properties, however, is not 
yet understood. One way to look at this question is by 
computational means using simulation techniques such 
as the discrete element method (DEM). This method 
takes into account the particulate nature of a colloid 
and allows for an investigation of the force distribu- 
tion inside the particle network during deformation as 
a function of the colloid's microstructure. However, this 
method needs to be supplied with initial particle config- 
urations. In preceding works, Brownian dynamics simu- 
lation (BD) was used to study the coagulation dynamics 
and the evolving microstructures in colloidal suspen- 
sions [6l[7]- These simulations were based on physical 
laws and widely accepted theories such as the Stoke's 
drag force, Brownian motion and the DLVO-theory [8], 
describing the inter-particle potential. The resulting mi- 
crostructures agree well with experiments [9] and can be 
used as initial particle configurations for further DEM 
simulations to establish the link between microstruc- 
ture and macroscopic mechanical properties. However, 
the BD method requires evaluating complex equations 
at each time-step in order to determine the various 
forces acting on the particles. Therefore, it is time- 
consuming, especially in the case of a repulsive energy 
barrier and moreover, structures with volume fractions 
exceeding 0.4 have not yet been simulated. For process- 
ing reasons, ceramic engineers are interested in prefer- 
ably high solid's phase volume fractions and in partic- 
ular in volume fractions exceeding 0.4. 

Inspired by the generation of heterogeneous micro- 
structures using ASP we developed the "void expan- 
sion method" (VEM), which allows for a fast and ef- 
ficient computational generation of porous microstruc- 
tures over a broad range of volume fractions and espe- 
cially those exceeding 0.4. 

In this publication, VEM is presented and the in- 
fluence of various simulation parameters such as the 
system size, the number of particles within the system 
or the elastic properties of the particles on the develop- 
ment of the mean coordination number is analyzed for 
a wide range of volume fractions between 0.2 and 0.55. 

2 Materials and Methods 

2.1 Discrete Element Method 

VEM is implemented using DEM [TU] and in particu- 
lar, the particle flow code in three dimensions (PFC 3D ) 



from Itasca Consulting Group, Inc., Minneapolis, Min- 
nesota, USA [TT] is used. DEM is an iterative method 
in which discrete spherical particles are used to build 
up more complex structures. At each point in time the 
forces on each particle are calculated. The time-step is 
chosen small enough to assume a constant force during 
the time-step, which allows for the linearization of the 
equations of motion enabling an efficient calculation of 
the particles' next positions and velocities. 

The forces on the particles included in our model 
arise from a linear elastic contact law between the parti- 
cles and damping. In particular, no other forces such as 
long range forces between particles or gravity are con- 
sidered. PFC 3D uses a soft-contact approach, wherein 
rigid particles are allowed to overlap at contact points. 
The contact law relates the forces acting on two con- 
tacting particles, in our case, linearly to the relative 
displacement between these particles. The magnitude 
of the normal contact force F n is given by Eq. (Q]) 

F n = k n U n (1) 

where k n denotes the normal stiffness and U n the 
overlap. The shear stiffness k s relates an incremental 
displacement in shear direction AU S to the shear con- 
tact force AF S via Eq. ©. 

AF S = k s AUs (2) 

The linear elastic contact law is thus parameterized 
by its normal and its shear particle stiffness. 

Energy dissipation is introduced via a local damp- 
ing term similar to that described by Cundall [T^] • The 
damping force, characterized by its damping coefficient 
d, is added to the equations of motion and is propor- 
tional to the force acting on the particle. Thereby, only 
accelerating motion is damped and the direction of the 
damping force is opposed to the particle's velocity [TT] , 

Thus, the forces in our model are characterized 
by three microscopic parameters: the particle's normal 
stiffness, its shear stiffness and the damping coefficient. 
In this work the inter-particle friction coefficient was 
set to zero in order to allow the maximum particle re- 
arrangement during the void expansion. 

2.2 Void Expansion Method 

VEM relies on two distinct kinds of particles: "struc- 
tural particles" that constitute the final microstructure 
and "void-particles" that are only used during the gen- 
eration of the structure. For clarity purposes, the first 
ones will be referred to as structural particles or simply 
particles and the latter ones will explicitly be termed 
void-particles throughout this publication. The phys- 
ically relevant macroscopic parameters characterizing 
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the final microstructures are N$, the number of struc- 
tural particles, rg, their radius and <Ps, the volume frac- 
tion of the structural particles. Using these parameters 
the edge length I of the cubic simulation box with pe- 
riodic boundary conditions is calculated using Eq. (|3j). 



Table 1 Simulation parameters 



1/3 



(3) 



The Ns particles are randomly placed in the simu- 
lation box with an initial particle radius of rs/(m +1), 
thus (rn+ 1) times smaller than the final rs, with m be- 
ing the number of subsequent radius blow-up steps. We 
use m — 10 in our simulations. At each blow-up step, 
the initial particle radius is added to the current parti- 
cle radius, followed by an equilibration of the structure, 
until, after the m th step, the final particle radius r$ is 
reached. This cyclic growing of the particles is needed 
in order to achieve volume fractions higher than ap- 
proximately 0.35 without a considerable particle over- 
lap which represents high local stresses. 

In addition to the structural particles, Ny void- 
particles with an initial radius ry <C rs are randomly 
placed in the simulation box. We used ry = 0.005 rs- 
After the structural particles have reached their final 
size the radius of the void-particles is increased cycli- 
cally. At each cycle their initial radius is added to 
their current radius, thereby simulating the swelling of 
the ASP. After each incremental increase of the void- 
particle radius, relaxation steps are performed in order 
to equilibrate the microstructure. This iterative pro- 
cedure is repeated until the structural and the void- 
particles are densely packed and any further increase 
of the void-particle size leads to a compaction of the 
particles, which is reflected by an increase of the strain 
energy inside the microstructure. Before each increase 
of the void particle radius the positions of the structural 
particles are stored, which allows for a subsequent anal- 
ysis of the microstructure as function of pore size, i.e. 
the void-particle's radius. 

In this study, the mean coordination number CN 
of the structural particles alone is used to characterize 
the evolving microstructures during the expansion of 
the void-particles. In particular, the coordination num- 
ber of a structural particle is given by the number of 
neighboring structural particles with a separation dis- 
tance smaller than d e = (1 + e)2rs, with e = 0.01. 

The density of bulk alumina was taken for the den- 
sity of the structural particles ps- The void-particle 
density py was set to a ten times smaller value py = 
ps/10 in order to reduce the inertia of the void- 
particles. Table Q] compiles the simulation parameters 
used in this work. 



Parameter 


Symbol 




Number of particles 


N S 


4000, 8000 


Particle radius 


rs 


2.5 x 10" 7 m 


Normal structural particle stiffness 


k n,S 


10 2 , 10 3 N/m 


Shear structural particle stiffness 




10" 3 , 10" 2 N/m 


Number of void-particles 


N v 


400 - 16000 


Normal void-particle stiffness 


kn,V 


10" 5 - 10 3 N/m 


Shear void-particle stiffness 


k 3 ,v 


10" 9 - 10" 1 N/m 


Damping coefficient 


d 


0.7 


Volume fraction 


<f>g 


0.2 - 0.55 


Structural particle density 


PS 


3690 kg/m 3 


Void-particle density 


PV 


369 kg/m 3 
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Fig. 1 Mean coordination number CN as a function of the void- 
to structural particle radius ratio q using 8000 structural parti- 
cles, 2000 void-particles and a volume fraction of 0.4. 



3 Results and Discussion 

In this study, the mean coordination number CN is 
used to characterize the evolving microstructures dur- 
ing the expansion of the void-particles. The evolution of 
CN as function of the void- to structural particle radius 
ratio q — ry /rs is shown in Fig. 1 using 8000 structural 
particles, 2000 void-particles and a volume fraction of 
0.4. Also, the conventions of the nomenclatures used 
throughout this paper are shown in this figure. 

The curve presents the three distinct regimes typical 
to all curves analyzed throughout this study (Fig. []}: 
The initial stage (I) is characterized by a small slope 
and a low mean coordination number. The slope within 
the second stage (II) increases drastically and the curve 
shows an inflection point. In the third stage (III), the 
particles are densely packed and any further increase 
of the coordination number is due to the compaction 
of the particles reflected by a significant increase of the 
structure's intrinsic strain energy. The transition stage 
can be interpreted as a phase change between stage I, 
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in which the particles can move freely and stage III, in 
which the particles' movements are arrested. The inflec- 
tion point in the transition region (II) is used to char- 
acterize the various curves that have been simulated. It 
is defined by three parameters: the void- to structural 
particle radius ratio qi = rv,i/fs, with ry,i being the 
void-particle radius at the inflection point, the mean 
coordination number CN and the maximum slope, de- 
noted by ACN/ Aq\i. Further characteristic parameters 
for the various curves are the void- to structural parti- 
cle number ratio n = Ny/Ns and the targeted volume 
fraction <P$, which is the volume fraction of the struc- 
tural particles alone. 

In the following, sensitivity analyses show the influ- 
ence of various VEM simulation parameters, especially 
the void- and structural particle numbers (Sect. 13. ip 
and the targeted volume fraction (Sect. 13. 2[) . on CN. 
These simulations use a void-particle normal and shear 
stiffness of 10 2 N/m and 10~ 2 N/m, respectively and 
a structural particle normal and shear stiffness of 10 3 
N/m and 10 -2 N/m, respectively. In Sect. 13.31 the scal- 
ing behavior of CN above the inflection point as func- 
tion of the total volume fraction is analyzed. The in- 
fluence of the void-particle stiffness on the evolving mi- 
crostructures is investigated thereafter in Sect. 13.41 for 
two distinct structural particle normal stiffnesses: 10 3 
N/m and 10 2 N/m. The structural particles' normal to 
shear stiffness ratio was fixed at 10 5 . 



3.1 Influence of the void- and structural particle 
numbers 

Microstructures with a volume fraction 4>s = 0.4 have 
been generated using various numbers of structural and 
void-particles, Ns and Ny, respectively. Ny essentially 
regulates the void size. Indeed, the higher Ny is cho- 
sen, the fewer blow-up steps are necessary in order to 
densely pack the structural particles. The influence of 
Ny on the VEM was probed with Ny ranging between 
400 and 16000 for N s = 8000. The size dependency of 
our system was tested with additional simulations for 
N s = 4000 using 2000 and 4000 void-particles. 

The simulations show that the evolving microstruc- 
tures depend on the void- to structural particle number 
ratio n, but they are independent on the individual ab- 
solute numbers of Ny and Ns- Increasing n shifts the 
transition region (II) towards a smaller void- to struc- 
tural particle radius ratio q. Fig. [2] shows qi normalized 
by n" 1 / 3 as a function of n yielding a constant value 
given in Eq. 




Fig. 2 Void- to structural particle radius ratio (qi normal- 
ized by n" 1 / 3 ) at the inflection point (circles, left scale) and 
the corresponding differential increase in coordination number 
(ACN/Aq\i normalized by n -1 / 3 , triangles, right scale) as a 
function of the void- to structural particle number ratio (n = 
Ny/Ng). Open and filled symbols denote simulations with TVg 
= 4000 and 8000, respectively. 



The total volume fraction at the inflection point <Pt, i 
(void and structural particles) is calculated via Eq. |5]). 



--$ s {l + qtn) (5) 

Thus, a constant value of f^n 1 / 3 entails a constant 
<?T,i- Indeed, the mean total volume fraction at the in- 
flection point is 62.3 ± 0.9 vol% close to 64 vol%, the 
characteristic volume fraction of random close packings 
(RCP) [13]. 

Fig. [2] further shows the differential increase in co- 
ordination number at the inflection point ACN/Aq\i 
normalized by n -1 / 3 as a function of n yielding a con- 
stant value given in Eq. ([6]). 

ACN/Aq\,n 1/3 = 9.16 ±0.68 (6) 

The best power law fit of ACN/Aq\i as a function of 
n yields an exponent of 0.3 slightly below 1/3 used for 
normalization in Fig.[2j The fit is very good as indicated 
by the correlation coefficient R 2 = 0.98. 

In particular, the use of two distinct numbers of 
structural particles (Ns — 4000 and 8000) for two void- 
to structural particle number ratios (n = 0.5 and 1.0) 
shows that the evolution of the mean coordination num- 
ber is independent on the system size. Virtually identi- 
cal values for both qtri 1 ^ 3 and ACW/AgljTi 1 / 3 as func- 
tion of n were obtained as shown in Fig. [51 

Combining Eq. g]) and Eq. © gives Eq. © 



q,n 



1/3 



0.82 ±0.01 



(4) 



ACiVl, = 11.2 9 Ag|i 

which, after integration, results in Eq. 
CN = 5.6qf + K 



(7) 



(8) 
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Fig. 3 Mean coordination number CNi as function of the void- 
to structural particle radius ratio qi, both at the inflection point, 
for various void- to structural particle number ratios n. 



Fig. 4 Mean coordination number at the inflection point CNi as 
a function of the volume fraction of the structural particles 



with K a constant. Thus, CN is expected to scale 
as qf, which means that CN essentially scales with the 
surface of the void particles. CNi as a function of qi and 
its fit (dashed line) are shown in Fig. [3l The trend is 
reproduced, however, the quality of the fit is rather low 
(R 2 = 0.76). In particular, the two data points with qi 
close to unity present a noticeable deviation of the gen- 
eral trend of increasing CN for increasing . This case 
is remarkable as it corresponds to an approximately 
monodispersed binary mixture of the structural and 
void-particles. A more detailed investigation, especially 
in the region of qi ~ 1, would require much more simu- 
lation runs spanning a wider range of particle number 
ratios n, which goes beyond this publication's scope. 



3.2 Influence of the volume fraction 

For processing reasons of ceramic bodies via colloidal 
routes engineers are usually interested in volume frac- 
tions as high as possible. A minimum of 40 vol% was 
required to perform uniaxial compression test on struc- 
tures fabricated using the DCC process [2]. Further- 
more, the volume fraction is easily accessible experi- 
mentally and is therefore widely used as a comparative 
value for various experiments. Thus, the influence of 
the volume fraction of the structural particles was 
investigated by means of CN as a function of q for <Ps 
ranging from 0.2 to 0.55 and for n — 0.5 (N$ = 4000 and 
Ny = 2000). The total volume fraction at the inflection 
point <pT,i (void- and structural particles) is calculated 
using Eq. |5]). Again, <Px,i yields a constant value of 
61.5 ± 0.5 vol%, slightly below the RCP limit. Thus, 
<I>T,i constitutes an upper boundary for <Ps ■ Indeed, for 



<Ps approaching $T,i and n ^ 0, qi approaches zero and 
any expansion of the void-particles is prohibited. 

Fig. [4] shows that CN increases with increasing vol- 
ume fraction As shown above, <Px,i is constant for 
•Ps ranging from 0.2 to 0.55, which indicates that at 
the inflection point, the structural particles are equally 
dense packed for all <Ps- Hence, CN may be expected to 
be a constant value, which seems to contradict Fig. 2] 
Plausibility considerations based on geometry give a 
possible explanation for this ^-dependence of CN. 
The sum of @s and $y,i equals &T,i, and is constant. 
Hence, rising <P$ results in a lower <Py,i and vice versa. 
Because Ny is constant, <&y,i OIu y changes by the vari- 
ation of ry.i. An increasing <&y t i entails an increase in 
ry^i, and vice versa. Structural particles in contact with 
void-particles have a lower coordination number than 
those which are only surrounded by other structural 
particles because only contacts between structural par- 
ticles are considered by definition of CN. The num- 
ber of structural particles in contact with void-particles 
scales essentially with the total surface of the void- 
particles. The smaller the total void surface, the less 
contacts between structural and void-particles exist and 
the more structural particles are only surrounded by 
other structural particles. Hence, increasing <!>s results 
in increasing CN. 

Similar considerations were used by Kruyt and 
Rothcnburg [15] who found a linear dependence be- 
tween a particle's coordination number and its radius in 
the case of two-dimensional assemblies of polydisperse 
particles. Further studies are necessary to confirm a lin- 
ear dependence between CN and <P$ as suggested by 
the line in Fig. [4] 
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Fig. 5 Scaling behavior of the coordination number CN above 
the inflection point in dependence of the total volume fraction 
•Pt for various curves (open symbols) and corresponding power 
law fits (lines). 



Fig. 6 Mean coordination number CN as a function of the void- 
to structural particle radius ratio q for a void- to structural parti- 
cle number ratio n = 0.5 and for various void-particle stiffnesses. 
The arrows indicate the direction of increasing void-particle stiff- 
ness. 



3.3 Scaling behavior of CN 

In order to analyze the scaling behavior of CN above 
CNi (region III) as a function of the total volume frac- 
tion #t the various curves analyzed in Sect. 13.11 and 
Sect. 13.21 were fitted using a power law given in Eq. ([9]) 

{CN - CNi) oc (<Z> T - $T.,if (9) 

with the exponent (3 as fit parameter. A selection 
of these curves is presented in Fig. [5] for simulations as 
function of the void- to structural particle number ratio 
n (circles) and as function of the volume fraction of the 
structural particles <I>s (triangles). The corresponding 
power law fits are shown as lines. 

For the simulations with varying n, an exponent 
/3 n = 0.39 ± 0.04 was obtained. The simulations in 
dependence of @s yield /3^ s = 0.35 ± 0.03. For all 
fits a very high correlation coefficient R 2 > 0.999 was 
achieved indicating excellent fits. In particular, our sim- 
ulations suggest that the exponents f3 is independent of 
n and <Ps- The average over all simulations results in 
an exponent f3 = 0.37 ± 0.04. This exponent is below 
the value 0.5 found in literature [T6lfT7]. where however 
dense instead of porous structures were considered. 

3.4 Influence of the void-particle stiffness 

The void-particle stiffness essentially regulates the ex- 
tent of the overlap of a void-particle with structural 
particles or other void-particles. Indeed, for a constant 
compressive force, a lower void-particle stiffness allows 
for larger overlaps between a void-particle and its neigh- 
bors. Thus, the void-particle stiffness is an important 



parameter for the microstructure and its evolution. For 
this investigation, the void-particle normal to shear 
stiffness ratio k n y /k s y is kept constant at 10 4 , n was 
fixed at 0.5 and <Ps — 0.4. Fig. [6] presents the analysis of 
CN as function of q for various values of k n y ranging 
from 10~ 5 N/m to 10 3 N/m and for k n , s = 10 3 N/m. 
Additionally, this analysis was performed for k n ,s = 10 2 
N/m and k n y ranging from 10~ 5 N/m to 10 2 N/m. In 
the following, the results are presented as function of 
the dimcnsionlcss void- to structural particle normal 
stiffness ratio K n = k n y jk n ^. 

Two principal behaviors are observed for increas- 
ing void-particle stiffness and thus for increasing K n : 
firstly, the inflection point is continuously shifted to- 
wards larger void-particle sizes and secondly, the transi- 
tion between region I and III becomes "sharper" . These 
two observations will be elaborated in the following. 

The shift of the inflection point towards larger void- 
particle sizes for increasing K n is summarized in Fig. 
showing the continuous increase of qi for rising K n . In 
particular, virtually identical curves are obtained for 
the two distinct values of k n ^s- For small normal stiff- 
ness ratios up to K n = 10~ 4 the values for qi as function 
of K n follow a logarithmic law given in Eq. (fT0|) . 

q % = 0.02 In K n + 1.16 (10) 

The fit is very good as indicated by a correlation 
coefficient R 2 = 0.993. For K n > 10~ 4 , ft levels off 
at approximately qi — 1.022, which corresponds to a 
total volume fraction of 0.61, close to the RCP volume 
fraction value. 

The characteristic void- to structural particle nor- 
mal stiffness ratio at which the transition from the log- 
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arithmic law to a constant value occurs is given by the 
intersection point of the respective fits. This value is 
termed and is found at 5.4 x 10~ 4 . 

A sharper, i.e. more "step" -like, transition between 
region I and III is observed in Fig.[S]for increasing k n y. 
Mathematically, the "sharpness" of the transition ex- 
presses in an increasing slope of the curves at their in- 
flection point, i.e. a higher value of ACN/ Aq\i, as sum- 
marized in Fig. [5] as function of K n . As in the case of 
qi(K n ) identical curves are obtained for the two values 
of fc n ,s- The values of AC N / Aq\i can be very well fitted 
against the void- to structural particle normal stiffness 
ratio K n up to 10~ 4 using a power law (Eq. (fTTjO . 

ACN/Aq\i = 410.3iC 2 ( n ) 

The correlation coefficient is R 2 = 0.997. For -Re- 
values higher than 10 -4 the AC/V/Ag^-values level off. 

The small exponent of 0.2 might also suggest a loga- 
rithmic dependence between AC N / Aq\i and K n , how- 
ever, the logarithmic fit is of considerably lower quality 
(R 2 = 0.91) compared to the power law fit. 

The characteristic kink is found at = 3.5 x 10~ 4 
N/m, which approximately corresponds to the value 
found for qi(K n ). 

Eq. (fT0|) and Eq. (fTTj) describe the independence 
of qi and ACN/ Aq\i, respectively, for K n < 10~ 4 . 
These equations allow to express ACN/Aq\i as func- 
tion of qi yielding a relation of the form ACN/ Aq\i ~ 
exp(qi). Indeed, an exponential fit gives a very high 
correlation coefficient (R 2 = 0.992). The integration of 
ACN/ Aq\i(qi) allows to predict the function CNi(qi), 
which also results in an exponential function. The fit of 
the simulated data using this function is good as well, 
achieving a correlation coefficient of R 2 = 0.93. 



4 Summary and Conclusions 

In this paper we presented VEM which allows for an 
efficient and fast computational generation of porous 
microstructures using DEM. The development of VEM 
was inspired by the experimental generation of hetero- 
geneous colloidal microstructures using ASP. VEM is 
a stochastic method in opposition to earlier used BD 
simulation, in which the physical processes during the 
coagulation of the colloidal suspension were simulated. 
This numerical description of physical processes how- 
ever requires much computing time. Thus, from a com- 
putational viewpoint, VEM is less intensive than BD 
simulation as VEM only includes a linear elastic con- 
tact law and damping. 

VEM permits to investigate the evolving mi- 
crostructure as a function of the void particle size. In 
order to characterize the microstructure we used the 
mean coordination number. For all simulation param- 
eters used throughout this research, the mean coordi- 
nation number as function of the void- to structural 
particle radius ratio exhibits the same characteristic 
"step"-shape: a small slope with low mean coordina- 
tion number in an initial stage, a transition stage with 
a sharp increase of the slope and, in a final stage, a 
small slope with high mean coordination number. The 
transition can be seen as a phase change between an 
initial stage, in which the particles can move freely and 
a final stage, in which the particles' movements are ar- 
rested. In this final stage, the particles are jammed and 
any further swelling of the void-particles corresponds to 
an increase in the strain energy in the structure. The 
inflection point in the transition stage characterizes the 
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obtained simulation curves and enables a comparison 
between the various curves. 

The sensitivity study of the mean coordination num- 
ber by variation of the void-particle number, the struc- 
tural particle number, the targeted volume fraction and 
the stiffness of the void-particles leads to the following 
results: 

1. The variation of the void- and structural particle 
numbers and of the volume fraction of the structural 
particles reveals an inflection point at a constant 
total volume fraction of approximately 62%, slightly 
below the RCP limit. 

2. The total volume fraction of 62% constitutes an up- 
per boundary for <P$ and therefore for VEM. 

3. The mean coordination number as function of the 
void- to structural particle radius ratio depends on 
the void- to structural particle number ratio alone 
and was found independent of the system size. 

4. Structures with volume fractions ranging from 0.2 
to 0.55 were successfully simulated using VEM. In 
particular, volume fractions above 0.4 were reached, 
which is crucial for a further simulation of the uniax- 
ial compression of colloidal structures using DEM. 

5. Above the inflection point, the mean coordination 
number as function of the total volume fraction fol- 
lows a power law with exponent (3 — 0.37 ± 0.04. 

6. An increasing void- to structural particle stiffness 
ratio K n reveals a twofold influence on the evolu- 
tion of the mean coordination number as a func- 
tion of the void- to structure-particle radius ratio: 
firstly, the inflection point is shifted to higher void- 
to structural particle radius ratios and secondly, the 
slope at the inflection point is increased. 

7. For small K n , the particle radius ratio qi and the 
maximum slope ACN/ Aq\i at the inflection point 
are nicely fitted versus K n using a power and a log- 
arithmic law, respectively. For K n approaching 1 
the curves level off. The transition from a power 
and logarithmic law, respectively, to a constant 
value is found at a constant void- to structural 
particle normal stiffness ratio K% of approximately 
4.5 x 10~ 4 . In particular, the curves qi(K n ) and 
ACN/Aq\i(K n ) do not depend on the structural 
particle stiffness. 

The computational time needed to generate a VEM 
structure is determined by the number of void radius 
blow-up steps that are necessary to densely pack the 
structural particles. The number of void-particle blow- 
up steps essentially depends on n and <Ps and is in- 
creasing for decreasing n or 4>s- Thus, smaller val- 
ues of n or <&s result in larger computational times. 
The simulations have further shown that for decreas- 



ing void-particle stiffness, less void blow-up steps are 
needed in order to reach the inflection point. However, 
for k n y < O.OOlfc^s lower packing density at the inflec- 
tion point are obtained for decreasing void-particle stiff- 
ness. In PFC 3D , the time-step essentially depends on 
the particle's mass and its stiffness as A/mp/fc, where 
mp is the smallest particle mass and k the largest par- 
ticle stiffness in the system. Hence, for given particle 
densities the time-step is determined by the structural 
particle normal stiffness k n ^s as long < k 

For k n y > k n ,s, the time-step decreases resulting in 
longer simulation times. 

To summarize, the VEM allows for an efficient com- 
putational generation of porous microstructures over a 
wide range of volume fractions and the various relations 
analyzed in this paper predict the influence of the VEM 
simulation parameters on the microstructures. Further 
research comprises the influence of inter-particle fric- 
tion, which in this study was set to zero in order to 
facilitate at most any particle rearrangements during 
the expansion of the void-particles. 

We use VEM for an efficient generation of porous 
colloidal microstructures over a wide range of volume 
fractions for subsequent simulation of the mechanical 
properties using DEM. Towards this goal, the VEM- 
microstructures have to withstand the comparison with 
experimentally determined microstructures using for 
example confocal laser microscopy |18j or with struc- 
tures obtained by other simulation techniques accu- 
rately describing the physical processes during the coag- 
ulation such as BD simulations [6j[T] . In this respect, an 
agreement in the mean coordination number is a neces- 
sary but not a sufficient condition. Additionally, further 
structural characterization methods, such as the pair 
correlation function [6] , the common neighbor distribu- 
tion function [19] or the recently introduced straight 
path distribution 20J must be considered in order to 
quantitatively compare the microstructures obtained by 
VEM to those obtained experimentally or using other 
computational techniques. 
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